function [S, R, X] = dpmm_gibbs_demo(K, n, lalpha, Kmax)
% A simple demo on dpmm_gibbs
%
%   [S, R, X] = dpmm_gibbs_demo(K, n, alpha);
%   [S, R, X] = dpmm_gibbs_demo(K, n, alpha, Kmax);
%
%       This function demonstrates the Gibbs sampling over DP mixture
%       model with 2D Gaussian component models.
%
%       Inputs:
%       - K:        the actual number of components
%                   It can also be input in form of [K, Kp], 
%                   where Kp is the number of prior component. (Kp <= K)
%       - n:        the number of observations from each component.
%       - lalpha:   the log-weight of sampling new particle
%       - Kmax:     the maximum number of estimated components
%
%       Outputs:
%       - S:    the consoliated samples
%       - R:    the raw samples
%       - X:    the data matrix
%

% Created by Dahua Lin, on Nov 16, 2010
%

%% prepare model

if isscalar(K)
    Kp = K;
else
    Kp = K(2);
    K = K(1);
end

if nargin < 4
    Kmax = inf;
end

disp('preparing model ...');

d = 2;

pri_sig = 8;
prior = gaussd.from_mp('s', zeros(d,1), pri_sig^2, 'ip');
mu = prior.sample(K);
comp_sig = 1;

gm = gaussgm_gp(prior, comp_sig^2);

sig_t = 0.2;
inherits = gaussgm_gp(gaussd.from_mp('s', mu(:, 1:Kp), sig_t^2, 'ip'), comp_sig^2); 

priw = constmat(1, Kp, n);
priq = constmat(1, Kp, 0.8);

%% synthesize data

disp('synthesizing data ...');

X = zeros(d, n*K);
for k = 1 : K
    X(:, (k-1)*n + (1:n)) = gm.sample(mu(:,k), n);
end

%% do sampling

disp('Doing sampling ...');

R = dpmm_gibbs(gm, X, lalpha, 100, ...
    'Inherits', inherits, 'PriW', priw, 'PriQ', priq, ...
    'BurnIn', 800, 'Interval', 5, 'Kmax', Kmax, 'Display', 'sample');


%% consoliate samples

disp('Consolidating samples ...');
S = consolidate_mmsample(R, 0.05 * n);

for k = 1 : numel(S.Id)    
    theta = S.Theta(:, k);
    fprintf('[%d]: (w = %.4g): [%s]\n', ...
        S.Id(k), S.W(k), num2str(theta', ' %.4f'));    
end

disp(' ');



%% visualize

% plot data

figure;
plot(X(1,:), X(2,:), 'b.', 'MarkerSize', 3);

% plot result

hold on;
G = gaussd.from_mp('s', S.Theta, comp_sig^2);
plot_ellipse(G, 1, 'g-', 'LineWidth', 2);
plot_ellipse(G, 3, 'g-');

axis equal;

